Viral Evolution and Adaptation as a Multivariate Branching Process 
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In the present paper we analyze the problem of adaptation and evolution of RNA virus popu- 
lations, by defining the basic stochastic model as a multivariate branching process. The defined 
stochastic process turns out to be well suited to describe several aspects of RNA viral populations. 
We show that in the absence of beneficial forces the model is exactly solvable. As a result it is 
possible to prove several key results directly related to known typical properties of these systems. 
Moreover, new insights on the dynamics of evolving virus populations can be foreseen. 
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I. INTRODUCTION 

RNA viruses exhibit a pronounced genetic diversity [l| . 
This variability allows RNA virus to better adapt to envi- 
ronmental challenges as represented by host and therapy 
pressures 0]. Due to the lack of a proofreading activity 
of viral RNA polymerases (average error incorporation 
rate in the order of 10~^ per nucleotide, per replication 
cycle Q), short generation times and huge population 
numbers, RNA viral populations may be viewed as a col- 
lection of particles bearing mutant genomes. As a con- 
sequence of high mutation rates, frequencies of mutants 
depend not only on their level of adaptation but on the 
probability of being faithfully replicated during viral ge- 
nomic RNA synthesis. Several studies have looked at 
viral diversification processes as a contributing cause of 
disease progression and of therapeutic strategies short- 
comings including vaccine trials 0j B- It has become 
important to understand the process by which virus ac- 
quire diversity and the dynamics and fluctuations of this 
diversity in time. However, understanding viral evolution 
in vivo has proven to be a very cumbersome accomplish- 
ment due to the so many variables present in the interplay 
between virus and their hosts. To name a few; the host 
defense pressures as the innate and "cognitive" immune 
responses, the use of antiviral drugs, the turnover rate 
of virus populations composed by viral replication and 
clearance, the elevated mutational rates of RNA virus, 
and the possible existence of structured viral reservoirs 
in infected patients. It is also important to take into ac- 
count the size of the viral innoculums at the moment of 
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infection, how frequent viral populations undergo bottle- 
neck passages within a host, how differently each infected 
individual may react to an incoming virus and finally how 
many viral variants are tightly associated with differen- 
tial biological capabilities. 

Traditionally, in an effort to make the viral evolution 
process more palpable, several groups have addressed this 
subject from different points of view. There is a sub- 
stantial amount of publications that studied virus pop- 
ulations during their evolution in experimental settings, 
for instance, cell cultures @ , by challenging the virus 
with population bottlenecks [61, |7[, or the introduction 
of antiviral drugs Q, including mutagens, or another 
competing viral population. Experiment outcomes were 
evaluated using viral replication kinetics, the intensity 
and quality of the observed mutational spectra and virus 
survival/extinction as final parameters. Several other 
groups have studied the process of viral evolution away 
from the bench but using mathematical and computa- 
tional tools [gl-fl^. These models arc quite tractable 
but there is always the risk of oversimplification. To 
escape from oversimplifying the interplay between virus 
and hosts a model needs to incorporate a few hard rules 
based on previous experimental data which has been gen- 
erated by the whole community of investigators address- 
ing viral evolution. Based on other groups experimen- 
tal data and previous mathematical models put forward 
by other investigators as the one presented by Lazaro et 
al. [l^ we sought to study a stochastic model for virus 
evolution that would be able to describe some general as- 
pects of RNA virus evolution. Here, RNA viral evolution 
is described by a multivariate branching process during 
which each round of replication is accompanied by the 
introduction of a single point mutatio n p er genome in 
the viral progeny. Drake and Holland Uj back in 1999 
have inferred, based on limited data, a central value for 
the RNA virus mutation rate per genome per replication 
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of /ig 0.76 and suggested the rate per round of cell in- 
fection of /Ltg ~ 1.5. In 2010, Sanjuan et al. (l6j revisiting 
this theme by reviewing a list of previous publications 
encountered RNA virus mutational rates in the order of 




w 4.64 for the bacteriophage Q/3 
) and /ig « 1.15 for hepatitis C virus 



10-* to 10- 
(Batschelet et al. 
(Cuevas et al. |18j). 

It has been demonstrated that virus populations may 
be reduced at the moment of infection, and only a few 
particles are able to start a new infection process in naive 
hosts [l^, [l^l. Abrupt reductions on RNA viral popu- 
lations known as population bottlenecks may eliminate 
population diversity and lead the virus to pathways to- 
wards extinction due to the exacerbated effects of genetic 
drift. An incoming virus population recovering from a 
transmission bottleneck event may show an asymptotic 
behavior resembling stationary equilibrium represented 
by the balance between two opposite forces classically 
identified with mutation and selection. This asymptotic 
behavior would occur if the environment is constant and 
enough time is allowed between two successive bottleneck 
events. The relaxation time between the bottleneck and 
the establishment of stationary equilibrium has been re- 
ferred to as the "recovering time" by Aguirre et al. [l3] . 

It has been pointed out (Drake and Holland [l^) that 
the basal value of RNA virus mutation rates is so large 
and RNA virus genomes are so informationally dense, 
that even a modest rate increase extinguishes the popu- 
lation. The frequent appearance of overlapping reading 
frames and multifunctional proteins augments the risk 
of a random mutation to have a deleterious impact and 
even more, multiply the effect of deleterious mutations. 
For example, the fraction of deleterious mutations out of 
random mutations occurring in vesicular stomatitis virus 
is around 70% (Sanjuan et al. [2l[). If the introduction 



of a mutagen to a replicating virus population is able 
to cause its extinction by increasing mutational rates, 
the process is known as chemical lethal mutagenesis and 
has been demonstrated in a number of viruses including 




23[ , human im- 
oliovirus type 
mphocytic 



the vesicular stomatitis virus (VSV) 
munodeficicncy virus type 1 (HIV-l) 
1 [23 , foot-and-mouth disease virus 
choriomeningitis virus p6j . Hanta virus |27j and Hep- 
atitis C virus H^. Accordingly, in the model, increases 
on mutational rates, and more specifically, on the dele- 
terious component of the mutational spectrum are able 
to push viral populations towards extinction. Our re- 
sults corroborate with the study from Bull, Sanjuan and 
Wilke [13] by showing that the sufficient condition for 
lethal mutagenesis involves mutational and ecological as- 
pects. Bull et al. arrived at a conjectural criteria for 
lethal mutagenesis by a heuristic and intuitive approach 
of possible general applicability. By applying the branch- 
ing process theory to the evolution of RNA viruses the 
lethal mutagenesis inequality proposed by Bull et al. [l2| 
is rigorously proven here. Furthermore, we describe four 
distinct regimes of RNA virus populations: transient 
regime, stationary equilibrium, extinction threshold, and 



extinction through lethal mutagenesis. 

The approach we adopt to the problem of virus adap- 
tation and evolution allows us to prove several key results 
directly related to known typical properties of these sys- 
tems and to get new insights on the dynamics of evolving 
virus populations. 

We note that in previous works [O, [11 [H-Hlj the 
properties of phenotypic models are discussed starting 
from a mean field linear model described by a mean ma- 
trix without reference to the underlying stochastic pro- 
cess modeling the microscopic dynamics of particle repli- 
cation. Here we deduce the matrix of first moments from 
a generating function of the stochastic process. Although 
the two matrices happen to coincide, it is important to 
stress that only from the generating function of the un- 
derlying stochastic process it is possible to fully discuss 
the validity of the model. 

Structure of the Paper. In section |TT] we describe a 
class of models for viral evolution and show that they 
define a multitype branching process. We explicitly com- 
pute the generating function and derive some elementary 
properties. In section lHll wc briefly recall some basic con- 
cepts and results from the theory of multitype branching 
processes. In section IIVI wc solve the spectral problem 
for the mean matrix of the model which allows us to ap- 
ply the results from section Hill Finally, in section IVl we 
present our conclusions and directions for future research. 



II. PHENOTYPIC MODELS FOR VIRAL 
EVOLUTION 

In this section wc describe a model for viral evolution 
and will show that it is naturally represented by a mul- 
tivariate branching stochastic process. Our description 
is along the same line of the probabilistic model intro- 
duced by Lazaro et al. [T^]. We interpret the notion of 
mutation probability as a general effect of probabilistic 
nature acting on the replication capability of individual 
viral particles, considered here as a measure of the par- 
ticle's fitness characterizing its phenotype. This effect is 
summarized by the definition of a stationary probabil- 
ity distribution which is used to set up a Galton- Watson 
branching process (Watson and Galton [s^) for the tem- 
poral evolution of the viral population. This probabil- 
ity distribution gives appropriate parameters to classify 
the asymptotic behavior of the viral population and to 
describe some of the non-equilibrium properties of the 
model. 

In other related publications the concept of mutation 
is extensively used as the cause of replication capacity 
change. Understanding that those changes constitute an 
observable output due to many different factors (of ge- 
netic and non-genetic nature), we prefer to use the gen- 
eral term "effect" over the replication capacity to char- 
acterize the three possible changes (deleterious, benefi- 
cial and neutral) that may happen with the viral particle 
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when it replicates. The precise definition of the three 
types of changes are given in the next section. 



A. Definition of the Model 

A number of viral infections starts with the transmis- 
sion of a relatively small number of viral particles from 
one organism to another one. The initial viral population 
starts replicating constrained by the unavoidable interac- 
tion with the host organism and evolves in time towards 
an eventual equilibrium. Each particle composing the 
population replicates in the cellular context that may dif- 
fer from cell to cell. Moreover each particle has different 
replication capabilities due to the natural genomic diver- 
sity found in viral populations in general. Therefore, it 
is reasonable to consider the viral population as a set of 
particles divided in groups of different replication capa- 
bilities measured in terms of the number of particles that 
one particle can produce. Each of those groups we call 
a class; the replication capability of a viral particle is an 
output of the process of interaction of that particle car- 
rying its genetic information with the cell environment. 
The replication capability is considered as a phenotypic 
character of the particle and therefore each class is con- 
sidered as a set of particles with a possible genotype di- 
versity expressing the same phenotype trait. The model 
we consider here does not take into account any informa- 
tion about the genomic diversity of any replicating class 
and therefore it should be classified as phenotypic model. 

We consider that the whole set of particles composing 
the viral population replicates at the same time in such a 
way that the evolution of the population is described as 
a succession of discrete viral generations. This assump- 
tion crucially depends on the clear definition of the time 
needed for a particle to replicate, referred by virologists 
as generation time. As it depends on the cell environment 
it is clear that this time period may vary from particle to 
particle replicating in different cells in such a way that 
the meaningful concept is a distribution of replication 
times with a possible clear mean value. The dispersion 
of the replication times can be considered small if we re- 
strict ourselves to homogeneous cell populations. Under 
these conditions we consider that no particle can be part 
of two successive generations. The possible impact of a 
subset of non replicating particles on the dynamics of the 
viral population is left to further studies. 

Suppose that we have a population of viruses that start 
evolving from an initial set of particles (population at 
t = 0), which is partitioned into classes according to the 
replication capacity of each particle, that is, where each 
particle of class produces no copies of itself, each parti- 
cle with class 1 produces one copy of itself, and so on. We 
assume that there is a maximum replication capacity R 
imposed by the natural limiting conditions under which 
any particle of the population replicates. Moreover, as 
the process of replication is controlled by chemical reac- 
tions involving specific enzymes and the template, it is 



reasonable to assume a mean bounded replication capac- 
ity per particle that is possibly typical for each specific 
virus. 

In the process of replication of a viral particle errors 
may occur at each replication cycle in the form of point 
mutations with possible impact on the replication ca- 
pacity of the progeny particles. Due to the intrinsic 
stochastic component of chemical reactions it is natu- 
ral to treat this point mutational cause as probabilistic. 
Another possible cause of change in the replication ca- 
pability in the viral offspring is clearly related to the 
cellular environment where the replication process takes 
place. As a result the time evolution of viral popula- 
tions should be viewed as a physical process strongly 
infiucnccd by stochasticity. Therefore we consider that 
the combined action of genetic and non-genetic causes 
may produce basically three types of replicative effects 
namely: deleterious with decreased replication capabil- 
ity from the parental to the progeny particles, benefi- 
cial when the replication capability increases and neu- 
tral with no change of the particles' replication capacity. 
In the model under consideration we define these three 
effects as: 

• deleterious ejfect: the replication capacity of the 
copied particle decreases by one. Note that when 
the particle has capacity of replication equal to it 
will not produce any copy of itself. 

• beneficial ejfect: the replication capacity of the 
copy increases by one. If the replication capacity is 
already the maximum allowed then the replication 
capacity of the copies will stay the same. 

• neutral ejfect: the replication capacity of the copies 
remain the same as the replication capacity of the 
original parental particle. 

For each type of effect we associate a probability at the 
particle scale applicable to every single replication event: 
d for the probability of the occurrence of a deleterious 
effect in one replication cycle of each particle, b for the 
probability of the occurrence of a beneficial effect in one 
replication cycle of each particle. The complementary 
probability c=l — 6 — dis the probability of the occur- 
rence of a neutral effect. In the case of in vitro experi- 
ments with homogeneous cell populations the parameters 
c, d and b may be considered as mutation probabilities. 

The simple phenotypic model is obtained by requiring 
that there are no beneficial effects in time, that is 6 = 0. 
This assumption is justified by several experimental re- 
sults. The frequencies between beneficial, deleterious and 
neutral mutations appearing in a replicating pop ulation 
have been already measured by prior studies |2ll . |331 - |40| . 
Taking their results together, it is reasonable to conclude 
that beneficial mutations could be as low as 1000 less fre- 
quent than either neutral or deleterious mutations. As a 
result the viral population would be submitted to a large 
number of successive deleterious and neutral changes and 
a comparatively small number of beneficial changes. 
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Prom what is described above it should become clear 
that the model assumes a scenario where a probabilis- 
tic processes at the cellular/ viral scale take place in the 
context of the interaction between the viral particle and 
the host cell. The combined effect of small scale pro- 
cesses are observed at the viral population scale in terms 
of collective (stable or not) properties. 

Based on the general aspects of the phenomenon of vi- 
ral replication it is compelling to to model it in terms of 
a branching process. In this perspective we define a dis- 
crete multitype Galton- Watson branching process for the 
evolution of the initial population, where the classes will 
be represented by the replication capabilities 0, 1, . . . , i?. 
The branching process is described by a sequence of 
vector- valued random variables {Zn ■ n g N} giving the 
number of particles in each replication class in the n-th 
generation. Thus Zn are vectors of non-negative integers 
satisfying the following assumption: if the size of the 
n-th generation is known, then the probability laws gov- 
erning the later generations does not depend on the sizes 
of generations preceding the n-th, that is the sequence 
{Zn : n S N} forms a markovian process. The initial 
population Z^ is represented by a vector of non-negative 
integers Zq = {Zq, Zq, . . . , Z^), which is non-zero and 
non-random. The temporal evolution of the population 
is obtained from a vector- valued discrete probability dis- 
tribution C = (Co, Cij • ■ • J Cfl) defined on the set of vectors 
with non-negative integer entries called the offspring dis- 
tribution of the branching process. For any vector with 
non-negative entries i ~ (z°, . . . , i^) one has that 

P(Z„+i=i|Z„ = e,) = Cr(i), (1) 

where = (0, . . . , 1, . . . , 0), with 1 in the r-th position. 
Thus, Cr{i) is the joint probability that an individual 
particle of class r (0 ^ r ^ i?) generates z° progeny 
particles in the class 0, progeny particles in the class 
1, . . . , i^ progeny particles in the class R. Note that 
any vector Z„ = {Z^, Zn, . . . , Zn) may be written as 
a sum J2r ^n^r and since each particle in Z„ may be 
seen as the initial condition of a new branching process 
independently of the others, equation ([IJ determines the 
probability laws for a general branching process as follows 

P(Z„+i ^ i\Zn = ErK^r) =l[Cr{if- . 

r 

In order to compute the offspring probability distri- 
bution C, for the simple phenotypic model, we start by 
observing that is non-zero only when i is of the form 
i = (0, . . . , i''"^, i'', . . . , 0) since a particle with replica- 
tion capability r can only produce progeny particles of 
the replication capability r or r — 1, moreover the en- 
tries i*""^ and i'' should satisfy i^^^ -\- i"^ = r. Thus we 
just need to compute the probabilities Cr on the vectors 
of the form ik = (0, . . . , r — fc, fc, . . . , 0). Suppose that a 
viral particle v with replication capacity r (0 ^ r ^ R) 
replicates itself producing new virus particles vi, . . . ,Vr. 
For each new particle v^, there are two possible outcomes 



regarding the type of change that may occur; neutral 
or deleterious, with probabilities c = 1 — d and d, re- 
spectively. Representing the result of the i-th replication 
event by a variable Xi, which can assume two values: if 
the effect is deleterious (failure) and 1 if the effect is neu- 
tral (success), the probability distribution of Xi is that 
of a Bernoulli trial with probability of occurrence of a 
neutral effect c = 1 — d (success), that is, 

P{Xi ^ k) ^ {1 - df d^-'' (fc = 0,l). 

The total number of neutral effects that occur when the 
original virus particle reproduces is a random variable 
Sr given by the sum of all variable Xi , since each copy is 
produced independently of the others, 

Sr — Xi + X2 + . . . + Xj- . 

That is, Sr counts the total number of neutral effects 
(successes) that occurred in the production of r virus 
particles vi,. . . ,Vr- It also represents the total number 
of particles that will have the same replication capacity 
r of the original particle v. It is well known (Feller [4l| ) 
that a sum of r independent and identically distributed 
Bernoulli random variables with probability c = 1 — 0? of 
success has a probability distribution given by the bino- 
mial distribution: 

P{Sr = fc) = binom(/c; r, 1 - d) = (^^ (1 - d)'' d""'^ . 

Since this is the probability that a class r virus particle 
V produces k progeny particles with the same replication 
capability as itself one has therefore 

Cr(0, . . . , r-fc. A:, . . . , 0) = P(5,. = fc) = binom(A:; r, l-d) . 

Given the offspring probability distribution C, one may 
set up a probability generating function f = (/o, . . . , Jb,) 
which is defined by the power series 

fr{zo, Zi, . . . , Zr) = ^ Crii) Zq . . . Zfl . 
i 

The probability generating function of the simple pheno- 
typic model is 

fo{zo,zi, . . . ,Zfl) = 1 

/i(zo, zi, . . . , zr) = dzo + czi 

f2{zQ,Zi,...,ZR)^{dzi+CZ2f ^2) 

/fl(^o, zi, . . . , zii) = [dzR^i + czr)^ 

Note that the functions are polynomials whose co- 
efficients are exactly binom(A:;r, 1 — d). This function 
completely determines the branching process. 

Now it is easy to obtain the general case where the 
beneficial effects have a non-zero contribution b. In this 
case, the binomial distribution is replaced by a trinomial 



5 



distribution (see Feller [41|) and the probability generat- 
ing function of the general phenotypic model is 

fo{zo,zi,...,ZR) = 1 

fi{zo, zi,...,zr) = dzo + czi + bz2 

f2{zo,zi,...,ZR) = {dzi +cz2 + bz3)^ 

: 

fn-iizo, zi, . . . , zr) = (dz_R-2 + C2fi-i + bzR)^^^ 
fnizo, zi,..., zr) = [dzR-i + (c + b)zR)^ 

Remark 1 It is worth to mention other variations of 
these models that share the same essential properties and 
are more adequate in different contexts. 

With Zero Class: In this variation, which is the ver- 
sion deduced above, particles of class ?- = are 
generated by the particles from class r = 1 . 

Without Zero Class: In this variation, the particle 
class is omitted and thus the probability gener- 
ating function has R variables and R components: 
omit the variable zq, the Rrst component fo and 
define /i(zi, . . . , zr) = + czi + bz2. Particles of 
class r = 1 undergoing a deleterious change are 
eliminated in the next generation. 



The evolution of the averages (Zn) of Z„ is given by 

=E(Z„|Zo) =M"Zo. (4) 

From the generating functions ^ and ^ it is trivial 
to compute the mean matrix of the phenotypic model. 
For the simple phenotypic model it is 
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Note that it is an upper triangular matrix. In the case 
of the general phenotypic model it is 
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B. Basic Properties of the Phenotypic Model 

We start by recalling that, when calculating probabili- 
ties and expectations, there is no loss of generality if one 
considers only initial populations consisting of just one 
particle of class r (0 ^ r ^ i?), since the general case can 
be decomposed as a sum of independent processes with 
this kind of initial population. All the relevant properties 
of the model can be deduced with this simplification. 

Wc shall introduce the notation Zq = 1 for the condi- 
tion Zq ~ Br, which is the initial population consisting of 
one particle of class r and zero particles of other classes. 
Thus P{Zi = ijZp = 1) = Ci (^)- A basic assumption 
in the theory of branching processes is that all the first 
moments are finite and that they are not all zero. Then 
one may consider the mean evolution matrix or the ma- 
trix of first moments M = {M^ } which describes how 
the averages of the sub-populations of particles in each 
replication class evolves in time: 

My =E(Z1|Z^ = 1), Vz,j = 0,...,i?. 

In terms of the probability generating function one has 

M.,.§|(1,1.^.,1). 

Denoting by /' the jacobian matrix of / one may write 



which is a tri-diagonal matrix. It is clear that d is related 
to upper part of the matrix, b to the lower part and c to 
the main diagonal. 

The mean matrix of the phenotypic model can be 
viewed as the adjacency matrix of a directed weighted 
graph where the nodes represent the particle classes ac- 
cording to their replication capacity and the arrows rep- 
resent the effect of decrease or increase of the replication 
capacity due to the replication process (see FIG.[T]). 

(a) 




12 3 



(&) 



12 3 




12 3 



FIG. 1. Graphs of mean matrices, (a) Simple phenotypic 
model, (b) General phenotypic model. The arrows are num- 
bered according to which there occurs a deleterious effect (d 
- dashed arrows) or a beneficial efi^ect (b - dotted arrows) or 
neutral effect (c - solid arrows). 



M = /'(!), where 1 = (1, 1, . . . , 1) . 
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III. RELEVANT RESULTS FROM THE 
THEORY OF BRANCHING PROCESSES 

In this section we collect a few definitions and results 
from the theory of branching process that will be neces- 
sary in our analysis of the phenotypic model. 

A. The Mean Matrix of a Branching Process 

Consider a multitype branching process with off- 
spring probability distribution and probability generat- 
ing function /. Suppose that ^ has all its first moments 
finite and not all zero. Then conditioning on the elemen- 
tary initial populations .Zq = 1 on may define the mean 
matrix M = {M^ } of the multitype branching process 
Zn by 

M,, =E(Zj|Zi5 = I) Vi,j = 0,...,i?. 

In general, a multitype Galton- Watson branching process 
can be classified into decomposable and indecomposable 
according to which its mean matrix is reducible or irre- 
ducible, respectively. A non- negative matrix M = {Mij} 
(0 ^ i,j ^ R) is called irreducible if for every pair of 
indices i and j, there exists a natural number n such 
that (M^) . . > and it is called reducible otherwise (see 
Gantmatcher [42[ ) . There is another characterization of 
irreducibility in terms of the graph of the matrix. 

The graph S{M) of M is defined to be the directed 
graph on R nodes {0, 1, . . . , i?}, each corresponding to a 
type of particle, in which there is a directed edge leading 
from node i to node j if and only if Mij 7^ 0. A graph 
g(M) is called path connected if for each pair of nodes 
(i, j) there is a sequence of directed edges leading from i 
to j. A matrix M is irreducible if and only if g(iW) is 
path connected (see Meyer Q). 

A multitype Galton- Watson branching process is called 
positively regular if its mean matrix M is primitive, that 
is, M" is positive for some positive integer n. In par- 
ticular, a positively regular branching process is inde- 
composable, since a primitive matrix is irreducible (see 
Gantmatcher or Meyer [i^l). Positive regularity is a 
standard assumption in the study of multitype branching 
processes, as it opens up the way to apply the powerful 
Perron-Frobenius theory (see Harris [4J] or Athreya and 
Ney [H). 

Example 1 Tiie classification of the phenotypic model 
according to the irreducibility or reducibility of its mean 
matrix is the following: 

(i) In the version '^with zero class" the mean matrix 
(H]) or ([5]) will have the first colmnn hlled with ze- 
ros, that is, they are not primitive matrices and thus 
the corresponding branching processes are not pos- 
itively regular. Moreover, a quick look at the graph 
9{M) in FIG. Q] (bj shows that the process is decom- 
posable since the node corresponding to particles of 



type does not have a direct arrow leading to other 
nodes. In the case of the simple phenotypic model, 
the corresponding graph S{M) is shown FIG. Q] (aj. 
Note that there are no dotted arrows since the prob- 
ability of benehcial effects is and so the graph 
is totally path disconnected, in other words, each 
"path component" of the graph consists of exactly 
one node. 

(a) In the version "without zero class" the mean ma- 
trix of both models can be obtained from ^ and 
([5]) by removing the first row and the Rrst column. 
Now the general phenotypic model becomes posi- 
tively regular, since the node corresponding to par- 
ticles of class 110 longer exists. The simple pheno- 
typic model still is decomposable, even without the 
node corresponding to particles of class 0. 

B. Malthusian Parameter and Extinction 
Probability 

Let q{M) denote the spectral radius of M , that is, if 
Ai, . . . , Aij are the eigenvalues of M then 

q{M) = max{|Ar|} . 

Since M is a non-negative matrix, it has at least one 
largest non-negative eigenvalue which coincides with its 
spectral radius (see Gantmatcher [i^ or Meyer j43j). 
When the largest eigenvalue is pos itive we shall call it, 
following Kimmel and Axelrod [461, the malthusian pa- 
rameter m of the branching process (see also Jagers et 
al. 

The malthusian parameter of a multitype Galton- 
Watson branching process plays the same role as the 
mean of the probability distribution of the offspring in 
a simple Galton- Watson process and its name is moti- 
vated by equation which implies that g{M^) = m", 
the average population size increases or decreases at a 
geometric rate, in accordance with the "Malthusian Law 
of Growth" . 

Finally, it follows from the theory of non-negative ma- 
trices that there is a left non-negative eigenvector v and 
a right non-negative eigenvector u corresponding to the 
eigenvalue m: 

M = TO and M u = mu , 

which can be normalized so that 

v^u = 1 and 1*m = I , (7) 

where is the transposed of the vector v. Moreover, 
when M is irreducible the left and right eigenvectors are 
positive (see Gantmatcher [42| or Meyer [43|). 

Let 7 — (70, . . . , 7i?) be the vector of extinction prob- 
abilities 

7r = P{Zn — for some n\ZQ — 1) , 
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the probability that the process eventuahy become ex- 
tinct given that initiahy there is exactly one particle of 
class r. In general, when the initial condition is given by 
a vector of non- negative integers Zq = {Z^, Zq, . . . , Z^) 
the extinction probability is 

R 

P(Z„ = for some n\Z^) = J| 7j^« . 

A basic result of the theory of branching processes is 
that the vector of extinction probabilities 7 is the solu- 
tion in [0, 1]^ with smallest components of the equation 

/(7)=7, (8) 

where / is the probability generating function. Observe 
that 1 is always a fixed point of /, that is, a solution 
of equation ([5]). Therefore, if there is no other solution 
of equation (|5]) in the unit cube [0,1]^ then the process 
always has probability 1 to become extinct. 

The main classification result in the indecomposable 
case, states that there are only three possible regimes 
(see Harris [13] or Athreya and Ney [isf): 

(i) If 771 > 1 then ^ 7 < 1 is the unique stable fixed 
point of / in the unit cube [0, 1]^ different than 1 
and the branching process is called super- critical. 
Therefore, with positive probability, the population 
will survive indefinitely. 

(ii) If m < 1 then 7 = 1 is the unique stable fixed 
point of / in the imit cube [0, 1]^ and the branch- 
ing process is called sub-critical. Therefore, with 
probability 1, the process will become extinct in fi- 
nite time. 

(iii) If m = 1 then 7 = 1 is the unique marginal fixed 
point of / in the unit cube [0, 1]^ and the branching 
process is called critical. Here, the expected time to 
extinction is infinite, despite the fact that extinction 
is bound to occur almost surely. 

Unfortunately this theorem does not cover all the in- 
teresting cases, one important example for us being the 
phenotypic model for viral evolution. Nevertheless, one 
of the earliest results about decomposable branching pro- 
cesses is the generalization of the classification, due to 
Sevastyanov (see Harris [3| and Jii^ina In the gen- 

eral decomposable case, there is a fourth alternative iden- 
tified by Sevastyanov [3, SI] and in order to formulate 
this condition wc need to introduce another important 
concept. A multitypc Galton- Watson branching process 
is called singular if its probability generating function 
is linear without constant term: f{z) = Mz. In this 
case, there is no branching since each particle produces 
exactly one particle that can be of any class and the 
process is equivalent to an ordinary finite Markov chain. 
More generally, a decomposable process may have singu- 
lar path components. Two nodes i and j are said to be 
in same path component if there is a sequence of directed 



edges leading from i to j and a sequence of directed edges 
leading from j to i. This procedure defines a partition 
of the set of nodes into equivalence classes, called path 
components of the graph ^[M). We say that a path 
component C of S(Af) is a singular path component, if 
any particle whose class is in C has probability 1 of pro- 
ducing, in the next generation, exactly one particle whose 
class is in C . Equivalently, the component functions of 
the probability generating function corresponding to the 
classes in a path component C are linear functions of the 
variables corresponding to the classes in the path com- 
ponent C. In other words, the "part" of the probability 
generating function corresponding to the classes in C is 
that of a singular branching process. The existence of 
singular components is obviously an obstruction to ex- 
tinction, for instance, in a decomposable singular process 
all path components arc singular. In fact, the result of 
Sevastyanov states that if there is at least one singular 
path component then the branching process never become 
extinct, no matter what is the value of the malthusian 
parameter. 

Example 2 The graph corresponding to the general 
phenotypic model (FIG. [B (b)) have two path compo- 
nents: {0} and {1, 2,3,..., R}. In the simple phenotypic 
model (FIG. [I] ('aJJ, the path components are exactly the 
sets containing one node, {0}, {!}, . . . , {R}. From the 
expressions of the generating functions ^ and Q it 
is clear that there are no singular path components in 
any of the models - simple or general, "with zero class" 
or "without zero class". Moreover, the general pheno- 
typic model "without zero class" is positively regular. 
Therefore, the phenotypic model displays only the three 
regimes determined by the malthusian parameter, which 
depends on the values of the parameters b, c, d and R. 

It is important to stress that the regime of a multitype 
branching process can not be read from the mean ma- 
trix alone (i.c, the malthusian parameter). Essentially 
this happens because of the existence of decomposable 
branching processes with singular components. 

Example 3 Consider the following generating functions: 

g{z, w) (1/2 + l/2z^ {dz + cwf) , 
h{z, w) = (z, {dz ~\~ cif)^) , 

where < c, d < 1 and c + d = 1. They have the same 
mean matrix given by M = ( J |^ ) malthusian 
parameter is m = max{l,2c}. It is easy to solve the 
fixed point equation ([8]) in both cases and compute the 
respective extinction probability vectors (7i,72)- for the 
function g we have that 71 = 1 and 72 = d'^/c^ if ^ 
d ^ i and 72 = 1 if | ^ c? ^ 1. For the function h we 
have that 71 = 72 = 0. Therefore, the branching process 
defined by g becomes extinct if and only if c ^ 1/2 while 
the branching process defined by h never becomes extinct 
irrespective of the value of the malthusian parameter! 
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C. Asymptotic Behaviour of Surviving Populations 

According to the "Mahhusian Law of Growth" it is ex- 
pected that a super-critical branching process will grow 
indefinitely at a geometric rate proportional to m" and 



we would like to write 



Wn , where W„ is 



random vector with a finite "asymptotic distribution of 
classes" when n ^ oo. The formalization of this heuris- 
tic argument is due to Kesten and Stigum (see [i^ [H^l 
for the case of indecomposable multitype branching pro- 
cesses and [5l[ for the case of a general decomposable 
multitype branching processes). 

Let us first recall the result in the indecomposable case 
(see Athreya and Ncy [1^). Consider a super-critical 
branching process with m > 1 and suppose that the vec- 
tor valued random variable ^ satisfies a technical condi- 
tion called K esten- Stigum "C log C " condition (see Lyons 
et al. (5^ and Olofsson j53l|). which is always satisfied 
in our case, since the probability distribution of the off- 
springs has finite support. It is natural to define the 
normalized random vector Wn = Zn/rn^. This normal- 
ized random vector has a limit when n — > cxd, that is, 
there exists a scalar random variable W ^ such that, 
with probability one, 

lim Wn = Wu, 

where u is the normalized right eigenvector correspond- 
ing to the malthusian parameter m and 



E(W|Zo) = v^Zo 



the 



where v is the left eigenvector corresponding to 
malthusian parameter m. 

An important step in the proof of Kesten-Stigum the- 
orem is the Kurtz [54| convergence of classes theorem: 



lim 



(almost surely). 



(9) 



Combi ning this with the Perron-Frobenius theorem (see 
Meyer 43|) one obtains 



lim 

n— f 00 



lim 



(Zn) 
\{Zn)\ 



(10) 



where the convergence in the first limit is in probabil- 
ity. The approach adopted in [2§-[3l| relates only to 
the second equality involving the limit of mean values of 
equation ()10p . By explicitly considering the microscopic 
model as a multivariate branching process the equality 
of the two limits in equation (jlO|) is guaranteed. This re- 
sult may be useful for computational simulations of the 
model, since one may compute the eigenvector u by sam- 
pling the population and taking averages. 

The meaning of the Kesten-Stigum theorem is that the 
total size of the population divided by m", converges to a 
random vector, but the relative proportions of the various 
"classes" approach fixed limits. Since we are assuming 



that the process is indecomposable the normalized right 
eigenvector u = (mq, . . . ,uji) is positive and is normal- 
ized so that '^^Ur — 1, therefore it defines a probability 
distribution on the set of classes {0, . . . , R}. It is called 
the asymptotic distribution of classes of the multitype 
branching process. 

In order to extend these results to the case where the 
branching process is decomposable one should employ the 
Frobenius normal form of the mean matrix M, which is 
reducible in this case (see Gantmatcher [13]). Kesten 
and Stigum [sij shows that it is possible, by rearranging 
the rows and columns, to rewrite the mean matrix in 
a block upper triangular form in such a way that the 
diagonal blocks are irreducible square matrices associated 
to components of the decomposable branching process. 
By a component of a decomposable branching process 
we mean a subset of classes such that their associated 
nodes in the graph S(iVf) forms a path component. Let 
{Cfc : ^ fc ^ N} be the set of components of S{M) 
ordered according to which Ck -< Ci if there is a sequence 
of directed edges leading from some i € Ck to some j S 
C/. Given two components Ck and C; define the sub- 
matrix 

M{k, I) = Mij with i e Ck ,j e Ci . 

Then, for each k, the square sub-matrix M{k) = M{k, k) 
is the irreducible mean matrix of the sub-process 

Z„(fc) = {Z; : ieCk}. 

Now the order of the components Ck allows us to rear- 
range the rows and columns of M in such a way that 



M = 



/M(0) M(0,1) M(0,2) 
M(l) M(l,2) 
M(2) 



V 











. .. M{0,N)\ 

... M{1,N) 

... M{2,N) 

M{N) ) 



(11) 



Therefore, the sub-process Zn{k) "receives input" from 
the sub- process Zn{l), with k < I, throughout the sub- 
matrix M{k,l). Note that if the sub-matrices M{k,l) 
arc all zero then the branching process splits as a sum of 
independent indecomposable branching processes. 

Example 4 The matrices (O and © of the simple and 
the general phenotypie models, respectively, already are 
in the normal form 

(i) In the simple phenotypie model we have that 
N = R, with one-dimensional diagonal sub-matrices 
M{k) = (A:(l — d)), with one-dimensional sub- 
matriees M{k,k -f 1) = ((fc + l)d) and one- 
dimensional sub-matrices M{k, l)^Qifl>k + l. 

(a) In the general phenotypie model we have N ~ I, 
with the Grst diagonal sub-matrix A4'(0) = (O) (or 
M{0) = (1) for the hrst variation of the model), 
the second diagonal sub-matrix M{1) = Mij, with 
i,j = 1,...,R and M(0, 1) = (dO---0). 
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Now observe that if = 1 with i G Ck then for I > k, 
the sub-process Zn{l) = for all n ^ 0. That is, the 
branching process behaves as if the sub-processes Zn{l) 
for all / > fc did not exist. Since each non-zero diago- 
nal sub-matrix M{1) is irreducible, it has a largest posi- 
tive eigenvalue m{l) and then we may define the effective 
malthusian parameter of the sequence of sub-processes 
(Z„(0),...,Z„(fc)) to be 



mc{k) 



max{m(/)| . 

l<k 



The simplest case is when all m{l) are simple eigenvalues 
of their respective sub-matrices M{1) - they are distinct 
amongst each other - this is exactly the case for matrices 
dl]) and ©. 

In Kestcn and Stigum (sij the result about the asymp- 
totic behaviour of irreducible super-critical branching 
process is generalized to the reducible case. The main 
theorem applied to the case where all m{l) are simple 
eigenvalues of their respective sub-matrices M{1) states 
that if the effective malthusian me(fc) > 1 and the 
"C log C" condition holds then for the normalized random 
vector Wn{k) = Zn/ {me{k))" there exists a scalar ran- 
dom variable W such that, with probability one, 

lini Wn{k) = Wu{k), 

where u{k) is the normalized right eigenvector corre- 
sponding to the effective malthusian parameter rrieik) 
and 

E(W^|Zo) = t;(fc)*Zo, 

where v{k) is the left eigenvector corresponding to 
the effective malthusian parameter me{k). Moreover, 
Kurtz's convergence of classes theorem ([9]) still holds. 
But one should note that the normalized right and 
left eigenvectors are not positive anymore. In fact, 
v{k) may have negative entries, but only those asso- 
ciated to the components Ci with I ^ k, for which 
Zq is zero. The right normalized eigenvector is of the 
form u{k) ^ (mq, . . . , Ur, 0, . . . , 0), where (mo,...,^^) 
is the non-negative right normalized eigenvector of the 
sub-matrix corresponding to the the sequence of sub- 
processes {Z„{0), . . . , Z„{k)), and so is a probability dis- 
tribution. 



D. Critical Behavior and Regime Transition 

The critical state separates the super-critical and the 
sub-critical regimes where the branching process has two 
distinct behaviors in time and thus characterizes the exis- 
tence of regime transition with genuine critical behavior. 
In fact, the decay of correlation functions described in the 
next section for the case of the simplest model clarifies 
this point. 

Although in a critical branching process Zn 0, al- 
most surely, when n — >■ oo, one still may obtain a mean- 
ingful asymptotic law by conditioning on non-extinction. 



See Mullikin [55[ and Joffe and Spitzer 56| for the in- 
decomposable case and Foster and Ney 53] for certain 
decomposable cases. 

In the indecomposable critical case Zn grows at a lin- 
ear rate prop ortional to n (see Harris |44l | or Athreya 
and Ney [45[), and so one should consider the normal- 
ized random vector Yn = Zn/n. If the second moments 
are finite and the branching process is non-singular, there 
is a scalar random variable Y ^ such that 



lim Yn 

n—^oo 



Yu given that Z„ ^ , 



where u is the normalized right eigenvector correspond- 
ing to the malthusian parameter m and with convergence 
only in distribution, which is weaker than the almost 
surely convergence in the super-critical case. 



IV. THE SIMPLE PHENOTYPIC MODEL 

For the simple phenotypic model it is easy to compute 
the eigenvalues of the mean matrix M: 



A,. = rc = r(l — d) r = 0, , 



In particular, the malthusian parameter is the largest 
positive eigenvalue 

TO = q{M) =Xr^Rc = R{1 - d) . (12) 

Therefore we have the following immediate result. 

Theorem 1 The simple phenotypic model has three dis- 
tinct regimes. 

(i) If R{1 — d) < 1 then the branching process is sub- 
critical. That is, with probability 1, the virus popu- 
lation becomes extinct in finite time. 

(ii) If R{\ — d) > 1 then the branching process is super- 
critical. That is, with positive probability, the virus 
population survives and grows indefinitely at an ex- 
ponential rate proportional to to" when n oo. 

(Hi) If R{\ — d) = 1 then the branching process is criti- 
cal. That is, with probability I, the virus population 
becomes extinct but this may take an infinite time 
to happen. 

Proof. This is a straightforward consequence of the main 
result about the classification of a multitype branching 
process and equation (|T2|) . □ 

Thus theorem [T] provides a partition of the parameter 
space {{d,R) : d 6 [0,1], i? G N} of the simple pheno- 
typic model into two regions (see FIG. [2]). 

It is also important, specially in order to describe the 
asymptotic behaviour in the super-critical case, to know 
the left eigenvectors v and right eigenvectors u corre- 
sponding to the eigenvalue Xr Let us write the left and 
right eigenvectors in components as 

V ^ {vq,vi,...,vr) and m = (mq, wi, . . . , u_r) 

and satisfying the normalization conditions ([T])- 
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FIG. 2. Graph of the function R = 1/(1 - d) (in blue). The 
region below this curve corresponds to the sub-critical param- 
eters (d, R) and the region above this curve corresponds to the 
super-critical parameters (d, R). The curve itself corresponds 
to the critical parameters {d,R). 



(i) In the version "with zero class" the left eigenvector 
V is given by 

.= (^(0,...,0,1) 

and the right eigenvector u is given by 

Uk = (^^^ (1 - d)'' d^~^ = binoni(fc; R,l - d) , 

with fc = 0, 1, . . . 

(ii) In the version "without zero class" there is no com- 
ponents vq and uq. The left eigenvector v is given 

by 

.^(^(0,...,0,1) 

and the right eigenvector u is given by 
1 fB> 



(1 - df d 



k jR-k 



1 - rf-R \ k, 

with k = 1, . . . , R. 

It is interesting to note that the simple phcnotypic model 
is a "completely solvable" branching process in the sense 
that we may explicitly solve the spectral problem for its 
mean matrix independently of the numerical values of the 
parameters. 

Next we turn to the computation of the extinction 
probabilities jr- In this case, it is necessary to solve 
a non-linear system of polynomial equations: 



20 = 1 

zi — dzQ + (1 — d)zi 
Z2 = {dzi + (1 - d)z2y 

zr = {dzR^i + (1 - d)ziiy 



(13) 



This may be done in a recursive way, since the equation 
for zo is already solved zq = 1 and the equation for Zk 
depends only on Zk and Zk-i- Thus we get for i? = 0, 1, 2: 



70 = 1 

71 = 1 

72 = < 



dV(l - d)2 for ^ d ^ i 
1 for i d 5$ 1 



When i? ^ 3 the formulas become very complicated and 
when i? ^ 5 the equation may not even be solvable by 
radicals, but in general one may write 



7r 



f{d) for s$ d 4 
1 for d =^ 1 



where dc = and f{d) is a strictly increasing smooth 
function on [0, 1[ satisfying: (i) /(O) = 0, (ii) f{dc) — 1, 
(iii) f{d) < 1 for ^ d < dc and (iv) lim^^i /(d) = -l-oo. 
This expression suggests that the surviving probabilities 
LJr = 1 — 7r can be interpreted as an order parame- 
ter associated to the occurrence of a phase transition 
when the deleterious probability d attains the critical 
point dc = ^-jr-, which marks the transition from super- 
criticality to sub-criticality, 

I g{d) for ^ d ^ dc 







for dc ^ d ^ 1 



where g{d) = 1 — /(d) and thus satisfies: (i) g{0) = 1, 
(ii) g{dc) = 0, (iii) 5(d) > for < d < dc and (iv) 
limd_).i g{d) ~ —00. Observe that for a fixed numeri- 
cal value of d, the system of equations ([T^ can be eas- 
ily solved by numerical approximation using Newton's 
method. For instance, in FIG. [3] we show the curves 
for the surviving probabilities as functions of d for 
r = 2,3,4,5,6,7. 

The result shows that, with respect to ujr, the model 
has a critical behavior in complete analogy to a second 
order phase transition (see FIG. [3]). Therefore, the criti- 
cal properties of the model can be characterized by means 
of relevant critical exponents. 

Finally, it is not difficult to see that for fixed d, the 
numbers 7^ satisfy 1 ^ 72 ^ 73 ^ ■ . . ^ 7_r and therefore 
the extinction probability for a general initial condition 
Zq = (Zq , . . . , Zq ) may be estimated far from the critical 
deleterious probability dc = (i? — 1)/-R by 



P{Zn ~ for some n\Zo) 



where |Z^| = 



Zq and near d, 



7f °l , (14) 
= (i?-l)/i?by 



P{Zn = for some n\Zo) 



(15) 



It has been demonstrated that large population pas- 
sages are able to increase the adaptability of virus popu- 
lations [1^. On the other hand, small population pas- 
sages represented by bottleneck events are capable to 
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FIG. 3. Curves for the surviving probability ciJr(rf) as function 
of d for r = 2, 3, 4, 5, 6, 7. 



increase the risk towards viral extinction. Among the 
aspects of abrupt population reductions are the exac- 
erbated effects of drift that coupled with the MuUer's 
hatchet principle [ssj may lead to the random and pro- 
gressive lost of the best adapted virus in a population. It 
also has been suggested that large virus populations bear- 
ing a significant phenotypic diversity are more adaptable 
to environment fluctuations and robust. It is correct to 
assume that large initial virus populations colonizing new 
hosts may show better survival probabilities than popula- 
tions recovering from bottlenecks. In this way the size of 
the viral innoculums may have an impact in the survival 
rates of different virus populations. 

From now on we shall split the analysis of the simple 
phenotypic model according to which it is sub-critical, 
super-critical or critical. 

It is important to note that the existence of a clear 
cut between regimes of survival and non survival popula- 
tions by means of a critical state is directly related to the 
problem of lethal mutagenesis for viral populations. In 
fact, proposition (i) in Theorem [J is precisely the Bull, 
Sanjuan and Wilke conjecture [l3] also represented in 
FIG.H 



A. The Sub-critical Regime: Lethal Mutagenesis 

The first consequence of theorem [1] is a proof, in the 
context of this model, of the conjecture of Bull, Sanjuan 
and Wilke for lethal mutagenesis [l^l . 



Corollary 2 In the simple phenotypic model, the virus 
population becomes extinct in finite time, with probability 
1, if the product of the neutral effect probability {1 — d) 
with the maximum replication capacity R is strictly less 
than one: 

(1 - < 1 . 

The main conclusion here is that the existence of lethal 
mutagenesis depends on "genetic components" (muta- 
tional rates) and other additional deleterious effects (host 
driven pressures intensifications), as well as on strict 
"ecological components" , namely, the maximum replica- 
tion capacity of the particles in the population and on 
the initial population size. As a result the viral popu- 
lation may reach extinction by increasing the number of 
deleterious mutations per replication cycle, by decreas- 
ing the value of R in the population or by a combination 
of the two mechanisms. The mutational strategy is the 
basis of treatments using mutagenic drugs (see Grotty et 
al. [1]) that induce errors in the generation process of 
new viral particles reducing their replication capacity. A 
straightforward consequence of corollary [5] is that a single 
particle showing the maximum replication capacity R is 
able to rescue a viral population driven to extinction by 
mutagenic drugs. If it is assumed that RNA virus pop- 
ulations correspond to a swarm of variants with distinct 
replication capacities, for a therapy to become effective 
it is important that it will eliminate the classes repre- 
sented by particles with highest replication capacities. 
As a conclusion the higher the replication capacity of the 
first particles infecting the organism the larger should 
be the number of deleterious mutations (or effects) and 
therefore the larger should be the drug concentration. 
This can be a clear limitation for treatments based on 
mutagenic drugs. 



B. The Super-critical Case: Relaxation and 
Equilibrium 

In the super-critical regime, the population grows at a 
geometric pace indefinitely. Nevertheless, there arc two 
distinct phases that occur during this growth: a transient 
phase ( "relaxation" or "recovery time" ) and a dynamical 
stationary phase. 



1. Relaxation towards equilibrium. 

An important question concerning the adaptation pro- 
cess of a viral population to the host environment is the 
typical time needed to achieve the equilibrium state. As 
the equilibrium is characterized by constant mean repli- 
cation capacity an obvious criteria to measure the time 
to achieve equilibrium would be by the vanishing varia- 
tion of this variable as used in other studies (Aguirre et 
al. [l3l)- Nevertheless, this method is clearly subjected 
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to the limitations of numerical accuracy with evident dis- 
advantages if one wants a sharp and universal criterion 
to differentiate populations from the point of view of how 
fast a population can be typically stabilized in a organ- 
ism. 

Viral populations are commonly submitted to transient 
regimes. As pointed out earlier the infection transmission 
process represents the passage of a small number of par- 
ticles from one organism to another in such a way that in 
this process the viral population is submitted to a sub- 
sequent bottle-neck effect during spreading of viruses in 
the host population. In order to approach the problem 
of relaxation after a bottleneck process in a more sound 
basis the natural quantity to be considered is the char- 
acteristic time derived from the decay of the mean auto- 
correlation function. The temporal correlation function 
C{n) is typically of the form cxp(— an) and the decay 
rate is given by the parameter a. The natural way to de- 
fine a characteristic time T to achieve equilibrium is by 
setting T = 1/a. In order to find the characteristic decay 
rates one should consider the recursive application of the 
mean matrix M on the initial population: Zq M^Zq- In 
fact, it is enough to consider the canonical initial popu- 
lation Zq = en = (0, 0, 1). By direct inspection it is 
easily verified that the decay of correlations is typically 
exponential and given by 

C(n) = exp ( - log(i?(l - d)) n) , 

where m = R{1 — d) is the malthusian parameter. The 
decay rate is therefore given by a = log(i?(l — d)). 

Among others, one possible application of this result 
relates to the very initial phase of the infection process. If 
we consider that during this phase the host immune sys- 
tem has not been yet stimulated against the virus, one 
can assume that the deleterious effects would be solely 
represented by the viral intrinsic mutation rates. There- 
fore, the largest the value of R, i.e., the largest the repli- 
cation capacity of the initial viral particle the fastest the 
progeny auto-correlation decays and reaches equilibrium 
stabilizing the viral population; intuitively the parameter 
R defines the degree of virulence of the infection during 
the early stage of the infective process. The increment 
of deleterious effects plays an opposite role on the decay 
rates. In fact, as it will be shown below the closest the 
parameter d is to its critical value c?c more time is needed 
to achieve equilibrium. 

2. The Dynamical Stationary State. 

When the simple phenotypic model is super-critical 
and is initialized with exactly one particle in the class 
r {Zq = 1) the effective malthusian parameter is mo = 
Xr = rc = r{l — d) with corresponding normalized right 
eigenvector u{r) = [uo{r), . . . , Ur{r), 0, . . . , O), where the 
components Uk{r), with /s = 0, . . . , r, are 

itfc(r) = binom(/c; r,l-d)^ T ) (1 - d)'^ . (16) 



Therefore, the simple phenotypic model has R — 1 dis- 
tinct asymptotic distributions of types of particles, de- 
scribing R — 1 distinct dynamical stationary states, char- 
acterized by their asymptotic distribution of classes given 
by (|16p (up to a random scalar perturbation) , each one of 
these being achieved when the branching process is ini- 
tialized with exactly one particle in the class r (Zq =1) 
for r = 2, . . . , i?, respectively. Note that when r — 0,1 
the process is always sub-critical. 

Theorem 3 // the simple phenotypic model is super- 
critical with malthusian parameter m ~ R{1 — d) and 
starts with at least one particle of class R then, in the 
long run, the relative number of particles in each class 
reaches a stable stationary dynamical state and is (up 
to a random scalar perturbation) distributed according 
to the Binomial Distribution: binom(fc;i?, I — d), where 
k ~ 0, . . . , R are the replication classes. 

Proof. This is consequence of Kesten-Stigum results 
about the asymptotic behaviour of super-critical mul- 
titype branching processes and the computation of the 
normalized right eigenvector associated to the malthu- 
sian parameter m ~ _R(1 — d) given by equation (|16p . □ 

From theorem [3] we immediately obtain: 

• The mean replication capacity is 

E(m) = i?(l - d) . 

• The phenotypic diversity is 

Var(M) = Rd{l - d) . 

It is well accepted that the phenotypic diversity is an im- 
portant characteristic of the viral population intuitively 
related to the idea of population robustness [s^, [^O] ■ In 
fact, a homogeneous population would be less fiexible 
from the point of view of adaptation. The variance as- 
sociated with the stationary state can be understood as 
a natural quantity to measure diversity. It shows that 
the maximum value of the phenotypic diversity r/4 is 
reached if d = 1/2 for any value of r. If _R > 2 the varia- 
tion of the phenotypic diversity as a function of d shows 
that there are two different domains to be considered: 
below d = 1/2 the diversity is an increasing function of 
d. It implies that if the population has a typical value of 
d < 1/2 the effect of inducing an increment of d (for in- 
stance using mutagenic drugs) increases the phenotypic 
diversity. For 1/2 < d < dc this effect reverses and di- 
versity decreases with increasing d. This result raises 
the question if in normal conditions the viral population 
adapt to the host environment guided by a principle of 
maximum phenotypic diversity or if the environmental 
conditions simply contribute to fix one possible value of 
diversity for the population that may vary from one to 
another host organism. Interesting enough, the natu- 
ral deleterious mutations has been measured for certain 
viruses and, as shown in the TABLE HI they are close to 
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the value d = 1/2. In the first case one could preview 
that the set point of the viral disease should be invari- 
ant (or with small variation) for all hosts. On the other 
hand the second hypothesis leads to the idea of different 
responses to treatment depending on the initial value of d 
before the adoption of treatment strategies to improve d. 
At the present the two scenarios may apply to different 
type of viruses and this point clearly has to be decided 
experimentally. 



Virus 


Ud 


(1 - d) = e-^-i 


REF. 


VSV 


0.692 


0.500 




TEV 


0.773 


0.461 


[37] 


'l'X174 


0.72 - 0.77 


0.48 - 0.46 


m 


Q/3 


0.74 - 0.86 


0.47 - 0.42 


m 



TABLE I. Experimental results of deleterious mutation rates: 
(VSV) vesicular stomatitis virus, (TEV) Tobacco etch virus 
and ($X174, Q/?) bacterial viruses. 

Another important consequence of the above results 
concerns the efficiency of the use of mutagenic drugs. In 
the region d < l/(i? + 1) < 1/2 the viral population's 
most representative particle is the fittest one (class R). 
If we assume that the drug action is deeply influenced by 
drug transport coefhcicnts in different host tissues, it is 
important to be assured that local drug concentrations 
will still eliminate the set of class R particles. If d in- 
creases beyond 1/{R+ 1) the representative particle of 
the population is not anymore the fittest one but a set 
of particles from different replication classes. Therefore 
the main drug target represents a group of average repli- 
cating particles of a population with higher phenotypic 
diversity in which resistance drug mutants can be con- 
tained. In this case one would say that the viral popula- 
tion displays a kind of endogenous strategy to scape the 
deleterious action of the mutagenic drug. If we assume 
that deleterious effects are small in the early stage of 
the infection process we should expect that at this stage 
the drug efficiency would be maximum reinforcing the 
successful practice of post exposure therapy, currently 
adopted in the case of HIV infections [H^] ■ 



C. The Critical Case: Extinction Threshold 

The clearest way to characterize the time behavior of 
the viral population at or around the critical point is 
through the typical time T to approach equilibrium de- 
rived from the decay of correlations described above. The 
expression T = 1/ log(i?(l — d)) shows that at the critical 
point the equilibrium state is never reached, i.e., the de- 
cay to equilibrium is at least non-exponential. A scaling 
exponent characterizing the behavior of T in the neigh- 
borhood of the critical point dc can be easily obtained. 



The expansion around dc = {R — 1)/R gives 

r«(i-4)M-4r^ 

Although it is always possible to calculate intermediate 
distributions of progeny, it is quite easy to see that at 
the critical point the time evolution of densities never 
achieves an invariant density. Unlike in the super-critical 
regime, the relative number of particles in each class/sub- 
population is never stable. Nevertheless, our preliminary 
results concerning the dynamics of fluctuations show that 
the time variation of the numbers of particles in each 
separated class follows a pattern such that the variation 
observed in one class is rigorously the same observed in 
all the others. This indicate a high level of correlation 
between the classes in complete analogy with critical phe- 
nomena of many physical systems. We conjecture that 
in the critical regime the highly correlated classes in the 
population behave as an inseparable whole such that the 
notion of the population divided in separated classes be- 
comes meaningless. In other words the correlation be- 
tween classes makes them behave as if they constitute 
one unique class, which reminds one of the basic proper- 
ties of the error threshold in Eigcn's theory [63j . In fact, 
according to Eigcn, when mutational rates are increased 
beyond a threshold, infinite viral populations are not any- 
more able to retain its best adapted variants. At this 
critical mutation level, selection is overruled by mutation 
and all variants share the same fitness status. Moreover, 
populations at Eigen's error threshold do not become ex- 
tinct, but well defined replication classes cease to exit, 
as particles hazardously wander through the surface of a 
flat landscape. If in the super-critical case the notion of 
the mean replication capacity and therefore that of the 
"mean viral particle" exists defining a typical scale in the 
system, in the critical case this notion is absent. There- 
fore, in using branching processes to model the time be- 
havior of viral populations the concept of error threshold 
should be identified with that of criticality. In the same 
direction of reasoning, in terms of branching processes 
the existence of lethal mutagenesis should be identified 
with that of critical behavior of the model. It is worth 
to note that there is a correspondence between Eigen's 
model of molecular evolution and the equilibrium sta- 
tistical mechanics of an inhomogeneous Ising system (see 
Leuthausser [gJ)- For further discussions about the error 
threshold in connection to extinction see (see (65l. l66j ) . 

The critical behavior of the model can also be observed 
through the survival probability uj for d < dc as show in 
FIG. [3] Expansion of the survival probability around the 
critical point by means of functional equation (eq. (|13p ) 
gives directly 

R 

Lj ^ 2— \d — dc\ ■ 

dc 

It is interesting to note that the critical exponents of T 
and UJ are the same found in critical behavior of a large 
class of dynamical random networks (sec Riordan and 
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Warnke [67|. This fact is reminiscent from the deep re- 
lation existing between branching process and random 
network theory, where the survival probability function 
of a branching process is identified with the order param- 
eter associated to the emergence of the giant cluster in 
a dynamical random network. This fundamental obser- 
vation goes back to Karp [68| and more recently it has 
become the central technique in the study of more general 
models of random networks (69| . The relation between 
the two theories is certainly expected to bring important 
results in the future. 

Finally, it is noteworthy that here we talk about crit- 
icality of a process taking place in time, and therefore 
the term critical phenomenon (imported from equilib- 
rium statistical mechanics of space distributed systems) 
is used to highlight the fact that the survival probabil- 
ity behaves like an order parameter and the amount of 
deleterious effects quantified by d behaves control 
parameter that can be changed by external means. 

V. CONCLUSIONS AND OUTLOOK 

Using the previous theoretical model for virus evo- 
lution proposed by Lazaro et al. [l^ and Aguirre et 
al. [Tsl . [l^ as a starting point we show that virus evolu- 
tion can be described by an exact solvable multivariate 
branching process. By applying our approach we are able 
to identify crucial aspects of the dynamics of replicating 
viral populations on a sound theoretical basis. Among 
these several aspects we are able to demonstrate that - 
as long as the beneficial effects are close to zero ~ the 
two main driving features of a virus population are the 
maximum replication capacity and the fraction of the 
population not affected by deleterious effects. Based on 
this result we show that, as proposed by Bull et al., if 
the product between the above mentioned parameters 
m = R{1 — d) yields a value less than one the population 
undergoes extinction. On the other hand, if m = R{l — d) 
is greater than one and the environment is constant we 
show that the population will reach an asymptotic sta- 
tionary state characterized by the stability of the replica- 
tive classes. However, the time to reach the stationary 
equilibrium strictly depends on how intense is the delete- 
rious effect, more precisely, the higher d the longer is the 
transient phase and when d approaches its critical value 
dc the transient tends to infinity. 

According to our explicit formulas for the progeny dis- 
tribution, we demonstrate that virus populations maxi- 
mize their phenotypic diversity by replicating with d near 
1/2, for any value of R. We speculate that this might 
be a universal property for RNA viruses that replicate 
under high mutational rates. In this way by increasing 
their phenotypic diversity viruses augments their chances 
of survival escaping and adapting to environmental pres- 
sures. Maintenance of high mutation rates makes it dif- 
ficult for a population to retain their best replicative 



classes. As a consequence, the best adapted classes are 
not the most represented ones in the population, thus 
not characterizing a classical Darwinian evolution pro- 
cess. As far as branching process modeling of viral evolu- 
tion is concerned its critical behavior partially resembles 
the concept of error threshold in Eigen's theory of molec- 
ular quasispecies. In this regime the replicative classes 
lose their independence in the sense that they become so 
much correlated that the whole set of classes behaves as 
a single one. 

We also demonstrate that by keeping the deleterious 
effects constant the survival probability of a virus pop- 
ulation will depend on its initial population size. By 
increasing the population size at time zero we push the 
survival probability curves, in the region before the criti- 
cal point, towards one (see FIG.|3]and equations ([T4| and 
(Uni)). According to this result it can be speculated that 
virus with greater innoculums have a better chance of 
survival colonizing new hosts. Interestingly enough and 
in a frontal disagreement to the above observations it has 
been shown that only a limited number of particles, and 
in some cases even one particle, is enough to start a new 
infectious process in a host [l^, [2^ . However, according 
to the model and as discussed before, the R parameter 
determines the success of an incoming virus population 
because the corresponding value of dc is uniquely given 
by R. The present work suggests that minimum innocu- 
lums must have at least one particle with replicative ca- 
pacity large enough in order to survive in the new host. 
We speculate that those particles with maximum replica- 
tive capacity should constitute the effective innoculum 
described in Zwart et al. [l9j . In fact, the experimen- 
tal data about viral load in HIV early infected patients 
strongly suggests that the host deleterious effects over 
the viral population are minimal and increase after the 
onset of the immunological response [T^I- We note that 
the characteristic form of this data can be easily repro- 
duced by the model (see Castro [Tlf and manuscript in 
preparation [t^). 

Finally, it is important to mention that the close re- 
lation of the theory of branching processes (as used in 
the present work) and dynamical Erdos-Renyi type net- 
works indicates that the latter may be brought to bear in 
the modeling of virus populations. The relation between 
these two theories is undoubtedly a research avenue with 
promising potential to improve our knowledge of the dy- 
namical laws governing the evolution and adaptation of 
viral populations. 
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